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SOLUTION  OF  SYSTEMS  OF  COMPLEX  LINEAR  EQUATIONS  IN  THE  L*  NORM 
WITH  CONSTRAINTS  ON  THE  UNKNOWNS 

I.  Introduction 

The  numerical  solution  of  general  systems  of  complex  linear  equations 
in  the  L^,  or  Chebyshev,  norm  is  a  mathematical  problem  that  arises  in 
several  applications.  The  objective  is  to  find  complex  values  for  all  the 
unknowns  so  that  the  maximum  magnitude  residual  of  the  system  of  equations 
is  a  minimum.  The  unknowns  are  not  allowed  to  assume  any  complex  value 
whatever;  instead,  they  are  required  to  satisfy  convex  constraints  of  the 
form  that  can  occur  in  the  applications. 

To  be  specific,  let  n,  m,  and  A  be  integers  satisfying  ®  >_  I,  n  >_  I, 
and  A  >_  1.  Let  the  matrices  A  e  CnXin,  B  c  CnX'*',  and  the  row  vectors 
f  £  Cm  ,  g  c  CA  ,  a  e  Cn  ,  d  e  Rn  ,  and  c  e  R*  be  given.  Let  A^  and 
denote  the  j-th  columns  of  matrices  A  and  B,  respectively.  The 
problem  investigated  in  this  report  can  be  stated  as  a  non-linear  mathe¬ 
matical  programming  problem  in  the  following  manner. 


Problem 


subject  to: 


min  e 
_n 


eeR.zeC 


zAj  -  fj|  <  e  ,  j  «  1 . . 

zBj  "  |  <  Cj  ,  j  -  1 , . . .  ,A 

zj  -  «j|  <  dj*  i  "  — 


(1) 

(2) 

(3) 

(A) 


The  vector  of  unknowns,  z,  is  also  taken  to  be  a  row  vector.  Notational 
convenience  is  the  only  reason  for  using  row  instead  of  column  vectors. 
It  is  assumed  throughout  this  report  that  c^  >  0  and  d^  >  0  for  all 
indices  j.  If  any  c^  or  d^  were  equal  to  0,  the  corresponding 
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inequality  could  be  replaced  by  an  equation  involving  no  absolute  values. 
These  equations  could  then  be  used  to  produce  a  problem  of  the  same 
mathematical  form  as  (1)  -  (4),  but  having  fewer  variables.  Further  dis¬ 
cussion  of  equality  constraints  is  given  later  in  this  section. 

If  there  are  no  constraints  (3)  -  (4)  and  if  the  system  is  under¬ 
determined  (i.e.,  m  rank  A  £  n),  then  simpler  methods  than  the  one  pre¬ 
sented  in  this  report  suffice  to  find  all  the  solutions  of  the  problem. 
Adjoining  the  constraints  (3)  -  (4),  however,  can  give  a  theoretically 
interesting  problem  even  if  it  is  under-determined.  For  example,  if 
m  »  n  -  rank  A,  then  the  problem  (1)  -  (2)  has  a  unique  solution  which  can 
be  found  by  solving  the  system  zA  -  f.  But  if  this  solution  does  not  sat¬ 
isfy  the  constraints  (3)  -  (4),  then  it  cannot  be  the  solution  of  the  prob¬ 
lem  (1)  -  (4).  The  method  presented  in  this  report  can  be  applied  whether 
the  system  is  over-determined  or  under-determined  with  equal  ease. 

In  the  approach  investigated  here,  this  problem  is  replaced  by  a 
linearized  problem.  The  linearized  problem  is  a  linear  program  which  is 
generated  in  such  a  way  that  the  relative  error  between  its  (exact)  solu¬ 
tion  and  the  (exact)  solution  of  the  original  non-linear  problem  can  be 
estimated  without  knowing  the  solution  of  either.  Furthermore,  the  maxi¬ 
mum  relative  error  can  easily  be  made  as  small  as  desired  by  selecting  an 
appropriate  linearized  problem.  See  Theorem  2  below. 

The  starting  point  for  the  linearized  problem  is  the  following  simple 
fact.  Let  x  and  y  denote  real  numbers  and  let  1  denote  the  imaginary 
unit .  Then 

lx  +  iyi  ■  max  (x  cos  9  +  y  sin  9)  .  (5) 

I  I  0  <  9  <  2* 
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A  proof  using  the  Cauchy-Schwartz  inequality  is  trivial.  The  maximum  in 
(5)  is  now  approximated  by  using  a  finite  number  of  angles  in  the  interval 
(0,2x).  Let  p  >  3  be  an  integer,  and  let  the  set  0  •  {0p...,0  }  be 
given.  Then  we  define  the  linearized  absolute  value 

lx  +  iy |  »  max  (x  cos  0  +  y  sin  0).  (6) 

I  le  0c0 


The  linearized  version  of  the  problem  results  from  replacing  every  absolute 
value  in  (2)  -  (4)  with  linearized  absolute  values. 

Linearized  Problem  min  e  ( 7 ) 

n 

eeR.zeC 


subject  to:  zA.  -  f ,  a  <  £ 

3  j  o  ~ 


j  -  l,...,m  (8) 

j  -  1 . *  (9) 

j  -  1 . n  (10) 


The  same  set  0  is  used  in  the  inequalities  (8)  -  (10)  merely  for  the  sake 
of  convenience;  it  is  certainly  possible  to  use  a  different  0  set  for  each 
of  the  m+n+i  inequalities  (2)  -  (4).  Substituting  in  real  and  imaginary 
parts  for  all  complex  quantities  and  applying  the  definition  (6),  it  be¬ 
comes  clear  that  the  linearized  problem  is  a  linear  program  in  2n  +  1 
unknowns  with  (m+n+A)p  inequalities.  See  section  II.  Even  for  modest 
values  of  m,n,X,  and  p,  this  linear  program  becomes  very  large  quickly. 

For  example,  if  n  m  Jl  ■  20,  m  ■  100,  and  p  ■  64,  it  has  41  variables  and 
8960  inequalities.  Furthermore,  the  inequalities  (8)  and  (9)  are  100X 
dense  since,  in  the  applications,  the  matrices  A  and  B  are  normally  100X 
dense. 
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The  linearized  problem  always  has  more  inequalities  than  unknowns,  so 
we  solve  it  by  solving  its  dual  using  the  revised  simplex  method.  The 
special  structure  of  the  dual  of  the  linearized  problem  is  used  effective¬ 
ly  in  several  ways.  One  prominent  feature  of  our  algorithm  is  a  signifi¬ 
cant  reduction  in  computer  storage.  Instead  of  (m+n+l)(2n+l)p  storage  lo¬ 
cations  required  in  a  straightforward  application  of  the  revised  simplex 
method,  the  method  presented  here  requires  only  2n(m+Jt)  +  2p  storage  lo¬ 
cations.  Another  prominent  feature  of  the  algorithm  is  that  the  most  nega¬ 
tive  reduced  cost  coefficient  of  all  the  variables  in  the  dual  problem 
is  computed  with  only  0(mn  log^p)  real  multiplications,  instead  of 
the  O(mnp)  real  multiplications  normally  required.  A  third  prominent 
feature  of  the  algorithm  is  that  the  solution  for  large  values  of  p  is 
approached  through  smaller  values  of  p.  That  is,  we  first  solve  the 
linearized  problem  for  p  »  4,  and  the  solution  for  p  »  4  is  used  as  an 
advanced  start  for  the  linearized  problem  with  p  »  8.  We  continue  in  this 
fashion,  doubling  p  at  each  step,  until  a  specified  value  of  p  is  at¬ 
tained.  In  this  way,  problems  with  p  as  large  as  2048  have  been  success¬ 
fully  solved.  These  three  features  have  been  incorporated  into  a  FORTRAN 
program  which  seems  to  be  practical  for  modest  values  of  m,  n,  and  i. 
even  for  large  values  of  p. 

Equality  constraints  of  the  form  zH^  ■  e^  ,  j  ■  l,...,r  with 
1  <_  r  <_  n  can  be  added  to  the  original  problem  if  desired.  One  way  to 
handle  them  in  this  context  is  by  the  method  of  elimination.  If  these 
equations  have  rank  q  <  r,  then  some  q  of  these  unknowns  can  be  solved 


for  explicitly  in  terns  of  the  remaining  n  -  q  unknowns.  These  q  un¬ 
knowns  can  then  be  substituted  out  of  the  original  problem,  yielding  a  new 
problem  with  only  n  -  q  unknowns.  Moreover,  the  new  problem  has  the  same 
mathematical  form  as  the  original  problem,  so  it  may  be  solved  in  the  man¬ 
ner  proposed  here.  For  this  reason,  equality  constraints  are  not  included 
in  the  original  problem  statement. 

The  set  9  used  in  (6)  can  be  an  arbitrary  set  of  angles.  However, 
it  is  convenient  to  restrict  0  in  two  ways.  First,  the  number  of  points 
in  0  is  required  to  be  a  multiple  of  4  and  a  power  of  2;  that  is,  p  ■  2 
with  K  >_  2.  Second,  the  points  in  0  are  required  to  be  the  arguments  of 
the  complex  p-th  roots  of  unity;  that  is, 

0k  -  (k-l)2x/p,  k»l,2 . .  (11) 

These  restrictions  allow  the  doubling  of  p  as  mentioned  above  and  simul¬ 
taneously  facilitate  considerable  computational  efficiencies  in  determining 
the  incoming  basic  variable.  See  section  III  below.  Another  reason  to  re¬ 
quire  (11)  is  that  for  fixed  p  it  gives  the  tightest  upper  bound  in  the 
inequality  (12)  below.  Unless  specifically  stated  otherwise,  0  will  be 
assumed  to  satisfy  both  these  restrictions  throughout  this  report. 

There  is  one  hazard  in  replacing  the  original  problem  (1)  -  (4)  with 
the  linearized  problem  (7)  -  (10).  The  linearized  constraints  have  a 
larger  feasible  region  than  the  original  constraints,  so  it  la  possible 
that  the  linearized  problem  has  solutions  for  some  values  of  p  even  when 
the  original  problem  is  actually  infeasible.  The  feasible  region  of  the 
original  problem  is  approximated  more  and  more  closely  as  p  is  increased. 


so  Che  linearized  problem  must  ultimately  fail  to  have  a  solution  for  suf¬ 


ficiently  large  p  when  the  original  problem  is  infeasible.  If  the  orig¬ 
inal  problem  is  in  some  sense  "nearly"  feasible,  but  in  reality  is  in¬ 
feasible,  the  linearized  problem  may  possess  solutions  for  very  large 
values  of  p.  Thus  one  may  be  deceived  in  certain  problems.  An  alterna¬ 
tive  viewpoint  is  that  any  false  solution  obtained  in  this  manner  to  in¬ 
feasible  problems  actually  represents  a  "reasonable"  solution  to  a  poorly 
defined  problem.  Whether  or  not  this  view  is  sensible  depends  on  the  ap¬ 
plication.  An  example  is  given  later  in  section  V. 

It  is  worth  mentioning  that  the  constraints  (4)  are  merely  instances 
of  the  constraints  (3),  and  similarly  for  the  constraints  (10)  and  (9). 
Distinguishing  these  two  kinds  of  constraints  is  desirable  not  only  for 
the  expected  reasons  of  simplicity  of  form  and  economy  of  computer  storage, 
but  also  because  the  simpler  constraints  present  themselves  implicitly  in 
the  unconstrained  problem.  That  is,  if  the  constraints  (9)  and  (10)  of  the 
linearized  problem  are  eliminated,  and  the  dual  of  the  resulting  linear 
program  solved  by  the  simplex  method,  then  some  of  the  constraints  (10)  re¬ 
appear  in  the  guise  of  artificial  variables.  This  will  be  made  clear  later 
in  section  II. 

.  The  following  theorem  lists  some  basic  properties  of  the  linearized 
absolute  value.  The  inequalities  (12)  are  especially  useful. 


»•  ** 


Theorem  1.  Let  p  >  4  be  an  even  integer,  and  let  the  elements  of  the  set 
0  satisfy  (11).  Then,  for  u  e  C  and  v  c  C, 


i. 

Me  - 

0  and  | u | q  - 

■  0  if 

and  only 

if 

u  *■  0. 

ii. 

|  u  +  V  j 

|0  <  Me  +  M 

0* 

iii. 

Given 

a  complex,  j 

He  ’ 

| a |  •  |u| 

e 

for  all 

u  if  and  only 

if 

arg  u 

e  0  • 

iv. 

Me 

Mfl  .  Ii® 

1  *  0  a  +  0 

n 

Me 

■  0 ,  and 

1 

lim  1 

u  T  ♦  0  1 
n|0 

au  1  =0. 

I  n| 

V. 

M.s 

|u|  <_  |u|q  sec 

(n/p). 

(12) 

Proof.  Properties  (i)  to  (iv)  are  straightforward.  The  first  inequality 
in  (12)  follows  immediately  from  (5)  and  (6).  The  second  inequality  in 
(12)  is  easily  proved  by  visualizing  the  set  of  all  u  e  C  such  that 
|u|q  «  1  as  an  equilateral  polygon  of  p  sides  whose  inscribed  circle 
is  Juj  »  1.  This  concludes  the  proof.  (Note  that  if  the  hypothesis  is 
weakened  to  p  ^  3,  the  theorem  remains  valid  except  for  the  first  equa¬ 
tion  in  (iv).) 

Although  the  linearized  absolute  value  is  not  a  norm  on  C,  it  is  a 
quasi-norm  because  of  properties  (i),  (ii),  and  (iv).  See  [1,  pp.  30-32]. 
Also,  the  asymptotic  result 

2 

sec  —  -  1  +  — _  ,  p  -*■  ®,  (13) 

P  2p^ 

shows  that  to  attain  a  relative  accuracy  of  5  significant  digits  (i.e.,  a 
relative  error  less  than  .5  *  10  ^)  in  (12)  requires  that  p  1024. 

The  next  theorem  connects  the  solution  of  the  original  problem  with 
the  solution  of  the  linearized  problem. 
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Theorem  2.  Let  e"  e  R  and  z"  e  Cn  solve  problem  (1)  -  (A),  and  let 


e'  £  R  and  z'  c  Cn  solve  the  linearized  problem  (7)  -  (10).  Then 

e'  £  e“  <  e'  sec(x/p)  (14) 

and 

jz'Aj  -  f  j  £  e  ’  sec(it/p)  ,  j  *  l,...,m  (15) 

jz'Bj  -  g  j  <  Cj  sec(x/p)  ,  j  -  1,...,A  (16) 

|zj  “  aj |  <  dj  sec(u/p)  ,  j  =  1 , . . . ,n  (17) 


Proof .  Inequalities  (15)  -  (17)  are  all  established  in  the  same  way.  For 

example,  since  we  must  have  |z!  -  a . L  <  d, ,  it  follows  from  (12)  that 

I  J  J  ~  J 


|zj  ~  aj |  i.  |zj  ~  aj |q  sec(n/p)  <  dj  sec(x/p.) 
The  following  sequence  of  inequalities  establishes  (14): 


*  max 

1  z '  A , 

-  f  1 

1  i 

j  19 

£  max 

1  z"A . 

-  f,L 

1  j 

j  10 

<_  max 

1  z"A 

-  f  A  1 

!  j 

jl 

<_  max 

I  z '  A, 

-  f,| 

1  j 

jt 

<  max 

Iz’A. 

-  f  L 

1  i 

i  10 

sec(n/p) 

J  1 

■  e'  sec(u/p) 

where  the  max  in  all  cases  is  over  j  -  l,...,m.  This  concludes  the  proof. 
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This  theorem  proves  that  the  solution  of  the  linearized  problem  is  an 
approximate  solution  of  the  original  nonlinear  problem  (1)  -  (4).  Fur¬ 
thermore,  it  proves  that  the  maximum  relative  error  in  this  approximate  so¬ 
lution  can  be  made  as  small  as  desired  by  appropriate  initial  choice  of  p. 

The  original  problem  (1)  -  (4)  has  a  mathematically  straightforward 
solution  if  all  the  quantities  are  real-valued  instead  of  complex;  that  is, 
the  real-valued  problem  is  an  ordinary  linear  program  in  n+1  variables  with 
2(m+n+A)  inequality  constraints  which  can  be  solved  in  a  finite  number  of 
steps  using  extensions  of  existing  methods  (see,  e.g.,  [2],  [3]).  The 
complex-valued  problem  is  much  less  simple.  Eliminating  complex  arithmetic 
from  the  problem  by  substituting  in  the  real  and  imaginary  parts  of  all 
complex  quantities  yields  (after  squaring  all  the  constraints)  a  mathemat¬ 
ical  programming  problem  in  2n+l  variables,  with  a  linear  objective  func¬ 
tion  and  m+n+A  quadratic  constraints.  No  method  is  known  which  can  solve 
this  problem  in  a  finite  number  of  steps.  Since  it  is  a  convex  programming 
problem  and  all  functions  involved  have  easily  obtained  derivatives  of  all 
orders,  a  great  many  different  algorithms  are  potentially  applicable  for 
its  approximate  solution.  However,  the  only  reference  (4]  known  to  the 
author  which  explicitly  studies  the  problem  (1)  -  (4)  uses  a  feasible 
directions  method.  At  each  step,  a  linear  program  is  solved  to  determine 
the  steepest  feasible  descent  direction,  a  line  search  determines  the  step 
length,  and  special  precautions  are  taken  to  prevent  zig-zagging,  or  jam¬ 
ming.  Also,  a  convergence  proof  is  supplied.  This  method  is  not  pursued 
further  here. 
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The  problem  (i)  -  (4)  can  be  viewed  as  a  semi -infinite  program  (SIP), 
as  is  discussed  briefly  in  section  III.  The  SIP  formulation  of  the  problem 
consisting  of  the  objective  function  (1)  with  only  the  constraints  (2)  has 
been  studied  in  (51,  (6],  and  [7],  where  it  appears  in  the  form  of  complex 
function  approximation.  The  linearized  absolute  value  (6)  was  first  given 
in  (6).  Theorem  1  and  a  less  general  form  of  Theorem  2  were  first  given 
in  (71. 

II.  Solution  of  the  linearized  problem 

A  solution  algorithm  for  the  linearized  problem  for  fixed  p  is  dis¬ 
cussed  in  this  section.  Special  attention  is  directed  to  aspects  of  the 
problem  which  give  rise  to  interesting  special  structure  as  a  result  of  the 

underlying  complex  arithmetic  of  the  problem.  Several  useful  theoretical 

0 

results  are  interspersed. 

It  is  first  established  that  the  linearized  problem  (7)  -  (10)  is  a 

linear  program,  as  stated  in  the  Introduction.  Denote  the  real  and  imagin- 

R  I 

ary  parts  of  any  quantity  x  by  x  and  x  ,  respectively,  whether  x 
be  a  number,  a  row  or  column  vector,  or  a  matrix.  By  definition  (6), 

jzA  -  f  I  -  max  ((zA  -  f  )R  cos  0  +  (zA  -  f  )*  sin  0], 

1  J  J QeQ  J  J  J  j 

so  the  m  inequalities  (8)  evidently  are  equivalent  to  the  system  of  mp 

Inequalities 

(zAj  -  fj)R  cos  0  +  (zAj  -  fj)1  sin  0  ^  e  ,  0  e  9  ,  j  *  l,...,m  .  (18) 
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Since 


/  *  *  AR  RAR  M  rR 

( zA  -  f  j  )  -  z  Aj  -  z  Aj  -  f  j 

(zAj  -  fj)1  -  +  zXA*  -  fj  , 

it  is  convenient  to  write  (18)  in  the  form 


(19) 


.  R  I  , 
[z  z  e  ] 


AR  cos  9  +  A*  sin  9 
AR  sin  9  -  A1  cos  9 
-  1 


<  [fR  cos  9  +  f1  sin  9] ,  0  c 


(20) 


where  we  have  defined  1  e  Rm  to  be  the  row  vector  whose  components  all 
equal  one.  Thus  for  each  9  in  (20)  there  corresponds  m  inequalities. 
The  inequalities  (9)  and  (10)  may  be  treated  similarly,  so  the  linearized 
problem  is  a  linear  program  in  2n+l  variables  and  m+n+1  Inequalities  as 
claimed.  The  linear  program  can  be  written  explicitly  as  follows: 


Primal  Problem. 


min 

,  R  I  ,  „2n+l 

[z  z  e]  e  R 


,  R  I  , 
[z  z  e] 


(21) 
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subject  to:  £  >.  0  and,  for  each  9  £  0, 


,  R  I  , 
[z  z  e] 


AR  cos  0  +  A*  sin  9 
AR  sin  0  -  A*  cos  0 
-1 


BR  cos  9  +  B*  sin  9 
BR  sin  9  -  B*  cos  9 


I  cos  0 
I  sin  0 
0 


(22) 


<  [fRcos  0  +  f*sin  9 


R  I 

c  +  g  cos  9  +  g  sin  9 


R  I 

d  +  a  cos  9  +  a  sin  9] . 


nxn 

We  employ  the  notation  I  c  R  for  the  identity  matrix  and  to  denote 

either  a  zero  row  or  a  zero  column  of  length  k  1.  Whether  it  denotes  a 

row  or  column  will  be  clear  from  the  context. 

The  primal  problem  is  solved  by  solving  its  dual  using  the  revised 

simplex  method.  The  simplex  (Lagrange)  multipliers  for  an  optimal  basic 

# 

solution  of  the  dual  problem  provide  a  solution  of  the  primal,  assuming  it 
to  be  feasible.  The  dual  can  be  written  in  one  of  the  standard  linear  pro¬ 
gramming  formats  by  explicitly  adding  a  slack  variable,  denoted  Q,  which 
arises  naturally  in  this  problem. 


Dual  Problem. 


min 

S  £  Ria><P,  T  c  R*Xp 
W  e  Rn*P,  Q  c  R 


(fR 

cos 

9k 

+ 

f1 

sin 

V 

sk 

R 

1 

g 

cos 

0 

+ 

g 

sin 

e  ) 

T 

k 

k 

k 

R 

a 

cos 

9k 

+ 

I 

a 

sin 

V 

W 

k 

(23) 
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Applying  the  asymmetric  form  of  the  standard  duality  relationship  (see, 
e.g.,  [8,  p.  69]),  it  is  easy  to  verify  that  the  dual  of  the  dual  prob¬ 
lem  (23)  -  (24)  is  equivalent  to  the  primal  problem  (21)  -  (22).  An 
alternative  statement  of  the  dual  problem  is  given  at  the  end  of  this 
section. 

The  slack  variable  Q  plays  a  special  role,  as  seen  in  the  next 
result. 

Theorem  3.  Let  the  matrices  S’  ^  0,  W’  0,  T’  0,  and  the  real  number 
Q'  ^  0  denote  any  optimal  basic  feasible  solution  of  the  dual  problem 
(23)  -  (24).  If  Q'  >  0,  then  the  optimal  value  of  the  objective  function 
in  the  primal  problem  is  zero. 

Proof .  The  proof  uses  the  asymmetric  form  of  the  Complementary 
Slackness  Theorem  (see,  e.g.,  [8,  p.  77]).  Let  [z'^  z'*  e']  e  R^n+* 
denote  the  simplex  multipliers  corresponding  to  the  optimal  basic  solu¬ 
tion  S’,  W',  T',  Q'.  By  the  Complementary  Slackness  Theorem,  Q'  >  0 
implies  that 


13 


I 


[z'R  z' 1  e  * ] 


-  e*  »  0. 


“0“ 
n 

0 

n 

1 

This  completes  the  proof. 

Except  for  the  slack  variable  Q,  every  basic  variable  of  the  dual 
problem  is  uniquely  identified  by  specifying  the  matrix  to  which  it  belongs 
together  with  its  location  (row  and  column  number)  in  this  matrix.  The 
matrix  names  S,  T,  and  W  clearly  correspond  to  the  inequality  systems  (8), 
(9),  and  (10),  respectively,  of  the  linearized  problem.  The  row  number  of 
a  basic  variable  identifies  the  particular  constraint  which  gave  rise  to 
it.  For  example,  all  the  dual  variables  in  row  q  of  matrix  T  would  be 
eliminated  from  the  dual  problem  if  the  q-th  inequality  in  (9)  ware  de¬ 
leted  from  the  linearized  problem.  Similarly,  the  column. number  of  a  basic 
variable  identifies  the  angle  in  the  set  9  to  which  it  may  be  said  to 
correspond.  Thus,  adjoining  a  new  angle  to  the  set  0  causes  one  new 
column  to  be  augmented  to  each  of  the  matrices  S,  T,  and  W. 

The  revised  simplex  algorithm,  as  applied  to  the  dual  problem,  is 
defined  in  general  terms  as  follows: 

Step  1.  Determine  an  initial  basic  feasible  solution  of  the  dual  problem. 

Step  2.  Compute  the  simplex  multipliers  corresponding  to  the  current 
basic  feasible  solution. 

Step  3.  Determine  the  incoming  variable  by  selecting  the  variable  having 
the  most  negative  reduced  cost  coefficient;  terminate  if  all  re¬ 
duced  cost  coefficients  are  nonnegative  —  the  primal  problem  is 
solved  by  the  current  simplex  multipliers. 


14 


Step  4.  Compute  the  column  of  the  incoming  variable  in  terms  of  the 
current  basis. 

Step  5.  Determine  the  outgoing  basic  variable  by  a  ratio  test;  terminate 
if  the  dual  objective  function  is  unbounded  below  —  the  primal 
problem  is  infeasible. 

Step  6.  Update  the  basis  inverse  and  current  basic  feasible  solution  by 
pivoting,  and  return  to  Step  2. 


The  special  structure  of  the  dual  problem  has  its  strongest  influence  on 
Steps  1,  3,  and  4.  These  effects  are  outlined  next.  Discussion  of  other 


important  aspects  of  the  algorithm" are  postponed  to  section  IV. 

The  dual  problem  is  already  in  canonical  form  for  initiating  the 
second  phase  of  the  simplex  algorithm.  In  other  words.  Step  1  is  trivial. 
To  see  this,  we  need  only  show  that  a  (2n+l)  x  (2n+l)  identity  matrix  can 
be  assembled  from  the  columns  of  the  coefficient  matrix  of  (24).  One 
readily  available  column  is  the  column  corresponding  to  the  slack  variable 
Q.  The  remaining  2n  columns  correspond  to  dual  variables  which  are  the 
components  of  two  particular  W  columns.  Recalling  the  construction  of 


the  set  9,  we  have  0^-0  so  that  cos  9^  «  1  and  sin  0^-0.  Hence 
one  of  these  W  columns  can  be  taken  to  be  W^.  Similarly,  the  other  is 


Wl4p/4  Since  9l+p/4 


x/2.  Hence  the  initial  basic  feasible  solution  is 


15 


The  simplex  multipliers  corresponding  to  (25)  can  be  derived  in  the  stand¬ 
ard  manner;  however,  it  is  more  convenient  to  derive  them  in  a  special  way 
later  in  this  section. 

The  remark  was  made  in  the  Introduction  that  the  constraints  (10) 
would  be  present  implicitly  even  when  not  explicitly  included  in  the 
linearized  problem.  If  the  constraints  (10)  are  deleted,  then  all  of  the 
W  variables  must  be  deleted  from  the  dual  problem.  An  easily  recognized 
initial  basic  feasible  solution  would  not  then  be  available,  so  it  would  be 
necessary  to  incorporate  2n  artificial  variables  into  the  dual  problem. 
These  artificial  variables,  together  with  the  slack  variable  Q,  would 
constitute  an  initial  basic  feasible  solution.  It  is  easy  to  see  that  this 
solution  is,  in  fact,  indistinguishable  from  (25);  hence,  the  artificial 
variables  are  actually  the  two  columns  and  in  disguise! 

The  only  difference  is  that  the  artificial  variables  would,  of  course,  have 
zero  cost  coefficients. 

The  initial  basic  feasible  solution  (25)  is  highly  degenerate  because 
all  but  one  of  the  constant  terms  in  (24)  are  zero.  As  discussed  in  [9], 
it  is  in  problems  of  this  general  kind  that  cycling  in  the  simplex  al¬ 
gorithm  has  been  occasionally  observed  in  practice.  Such  cycling  was  ob¬ 
served  in  an  example  given  in  this  report.  However,  a  trivial  modification 
of  the  tie-breaking  rule  in  the  ratio  test  for  the  outgoing  basic  variable, 
together  with  "preferential  treatment”  of  certain  incoming  variables,  seems 
to  completely  avoid  the  problem.  Further  discussion  of  cycling  in  the  dual 
problem  is  postponed  to  section  IV. 
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The  cost  coefficients  and  the  columns  of  any  dual  variable  can  be 
found  by  inspection  of  (23)  -  (24).  It  is  better,  however,  to  think  of 
them  in  the  complex  arithmetic  format  given  in  Table  I,  which  is 


Table  1.  Dual  variable  cost  coefficients  and  columns  in  complex  arithmetic 
format. 

readily  verified.  Utilizing  this  format  makes  a  very  significant  reduction 
in  computer  storage  easy  to  understand.  Instead  of  storing  the  (m+n+i)p 
columns  of  all  the  dual  variables,  the  column  of  a  dual  variable  is  con¬ 
structed  only  if  it  is  needed.  As  indicated  in  Table  1,  this  requires 


storing  only  the  matrices  A  and  B  and  the  cosines  and  sines  of  each  of 
the  p  angles  in  0.  Assuming  the  necessary  2p  cosines  and  sines  are 
computed  once  and  for  all  at  the  outset,  only  n  complex  multiplications 
are  required  to  construct  the  column  of  any  specified  dual  variable. 

This  is  a  very  small  price  to  pay  for  reducing  the  storage  from 
( 2n+l)(m+n+i)p  words  to  only  2n(m+n+i)+2p  words.  Furthermore,  the  col¬ 
umns  of  the  dual  variables  W  ^  are  constructed  from  the  identity  matrix 
I,  which  need  not  be  explicitly  stored.  Hence  the  total  storage  necessary 
for  constructing  the  column  of  any  dual  variable  is  merely  2n(m+£)+2p 
words . 

A  very  efficient  method  of  computing  the  smallest  reduced  cost 
coefficient  in  Step  3  of  the  revised  simplex  algorithm  is  now  discussed. 
This  method  is  particularly  interesting  because  none  of  the  columns  6f  the 
dual  variables  are  explicitly  needed.  The  only  data  requited  are  the 
original  complex  matrices  A  and  B  and  the  sines  and  cosines  of  the 
angles  in  0  .  Let  X  be  any  real  row  vector  of  simplex  multipliers  for 
the  dual  problem;  thus,  X  is  of  length  2n+l.  The  vector  X  defines  a 
complex  row  vector  z  e  C°  and  a  real  number  e  e  R  by  the  identification 


,  A  ,  R 
X  ■  (z 


-e]  £  R 
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This  definition  is  reasonable  considering  Che  statement  (21)  -  (22)  of  the 

primal  problem.  The  reduced  cost  of  the  dual  variable  is  the  cost 

coefficient  of  S  .  minus  the  row  X.  times  the  column  of  S..  •  Using 
J*  jfc 

(26)  and  Table  1  shows  that 


„1k  A  "i0k  xR  r  R  I  , 

CS  *  (fje  )  “  L  2  2 


e  - 


(!Aj  -  f3>  • 


-i8. 


-19.  R 

(*,«  k  ) 

-10.  I 
-  ( A^e  k  ) 

1 


(27) 


so  the  minimum  reduced  cost  coefficient  of  all  p  variables  in  the  j-th 
row  of  S  must  be 


d  -  min  Cf 
l<k<p 


C  ->2Aj  •  fj 


»  j  ■  1»2,... ,m. 


(28) 


The  smallest  reduced  cost  coefficient  of  all  the  dual  variables  of  the 
matrix  S  is  then 

C,,  ■  min  C^,  -  c  -  max  IzA  -  f  I  .  (29) 

b  l<,j<m  &  K,j<m  1  3  J,0 


Similarly,  the  minimum  reduced  cost  coefficients  over  all  the  dual  vari¬ 
ables  of  T  and  W  are 
,  A 


min  (c  -  |*B4  -  gj  ) 


J 


(30) 


e 
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and 


s. 


A 


min 

Kj<n 


(31) 


respectively.  The  smallest  reduced  cost  of  all  the  variables  of  the  dual 
problem  is 


CSWT  “  111111  ^CS’  Sr  CT^' 


(32) 


The  actual  value  of  the  minimum  reduced  cost  CSWT  is  not  particular¬ 
ly  important.  What  is  Important  is  knowing  for  which  dual  variable  the 
minimum  is  attained.  This  requires  knowing  which  of  the  three  quantities 
Cg ,  Cy,  or  CT  is  smallest  as  well  as  knowing  for  which  index  j  the 
minimum  values  of  (29)  -  (31)  are  attained.  These  two  pieces  of  data  tell 
the  row  number  and  the  correct  matrix  name  of  the  incoming  dual  variable. 

To  determine  the  column  number,  we  must  also  know  which  angle  9^  c  0 
gives  the  largest  projection  (i.e.,  the  linearized  absolute  value  (6))  at 
the  minimal  index  j .  Since  the  angle  0^  may  not  be  unique  because  of 
possible  ties  in  (6),  a  consistent  tie-breaking  rule  must  be  defined.  This 
rule,  which  we  will  call  the  minimal  clockwise  index  (MCI)  rule,  is  import¬ 
ant  because  it  determines  unambiguously  the  name  of  the  incoming  dual  vari¬ 
able  in  the  simplex  method. 

Let  u  €  C,  and  let  Uq  be  the  set  of  those  angles  0  e  0  for  which 

the  maximum  in  (6)  is  attained.  There  are  three  cases.  First,  if  u. 

0 


has  precisely  one  element,  the  MCI  of  u  is  defined  to  be  the  index  of 
that  element.  Second,  if  u^  has  precisely  two  elements,  say  0^  and 
and  neither  k  or  j  equals  p,  then  the  MCI  of  u  is  defined  to  be 
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min  {k,j};  on  the  other  hand,  if  either  k  -  p  or  j  -  p,  then  the  MCI  of 
u  is  taken  to  be  p.  The  reason  for  the  exception  is  clear  from  the  geo¬ 
metry  of  the  problem.  Third,  if  has  more  than  two  elements,  then  it 

must  be  that  u  *  0  and  =  0 ,  so  the  MCI  of  u  is  defined  to  be  1. 


Thinking  of  the  simplex  multipliers  as  a  complex  vector  of  length  n 
and  a  real  number,  as  in  (26),  is  “natural"  in  this  problem.  The  complex 
number  zA^  -  f  ,  for  example,  can  be  calculated  entirely  in  complex 
arithmetic.  This  makes  any  computer  program  implementing  the  algorithm 
easier  to  write  and  also  probably  makes  it  execute  more  efficiently  (when 
coded  intrinsically  in  COMPLEX  mode).  Furthermore,  the  Idea  emerges  that 
the  computation  of  the  linearized  absolute  value  | zA^  -  f  and  its  MCI 
should  be  treated  as  a  separate  optimization  problem  having  its  own  special 


features.  This  subproblem  is  discussed  next. 

The  computation  of  the  linearized  absolute  value  and  corresponding 
MCI  must  be  undertaken  for  m+n+I  complex  numbers  during  each  iteration  of 
the  simplex  algorithm  in  Step  3.  A  brute  force  approach  using  the  defini¬ 
tion  (6)  might  require  as  many  as  2p  real  multiplications  for  each  com¬ 
plex  number.  Such  an  approach  is  quite  inefficient  and  does  not  exploit 
the  special  form  of  the  set  9.  For  p  -  4,  it  is  clear  that  comparison 
tests  alone  suffice  to  solve  this  subproblem.  For  p  ^  8,  we  claim  that 
comparison  tests  and  at  most  2  log^P  -  5  real  multiplications  are  suffi¬ 
cient.  To  see  this,  first  determine  the  octant  of  the  complex  plane  in 
which  the  given  number  lies,  and  then  determine  whether  it  lies  above  or 
below  the  45°  line  bisecting  the  octant.  This  can  be  done  using  comparison 
tests  only.  Now  that  the  "half-octant”  in  which  the  number  lies  is  known. 
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its  projections  onto  the  bounding  rays  of  this  half-octant  can  be  computed 
in  this  special  case  using  only  one  multiplication.  If  p  ■  8,  a  final 
comparison  test  ends  the  problem.  If  p  _>  16,  then  the  larger  of  the  two 
projections  reveals  the  "quarter-octant"  in  which  the  number  must  lie.  The 
projection  onto  one  of  the  bounding  rays  of  this  quarter-octant  is  already 
known,  so  it  is  only  necessary  to  compute  the  projection  onto  the  other 
bounding  ray.  This  requires  2  real  multiplications.  If  p  =  16,  a  final 
comparison  test  ends  the  problem.  If  p  >  32,  we  continue  as  before. 
Counting  the  total  possible  number  of  steps  proves  our  claim.  This  bi¬ 
section  method  works  because  we  have  required  the  set  0  to  contain  only 
the  points  (11).  Other  choices  of  0  would  require  different  methods. 

The  number  of  real  multiplications  required  to  complete  Step  3  of  the 
revised  simplex  algorithm  using  the  methods  discussed  above  is  significant¬ 
ly  less  than  that  required  in  the  usual  approach.  Given  the  simplex  multi¬ 
pliers,  the  straightforward  method  requires  the  computation  of 
(m+n+A)p  -  (2n+l)  real  inner  products  of  length  2n.  Taking  account  of  th«- 
simple  form  of  the  W  columns  gives  a  total  of  approximately 

(2p-4)[n(m+A)+l)  +  4n(m+A-n-l/ 2)  (33) 

real  multiplications.  The  special  methods  discussed  above  require  m+n+A 
complex  inner  products  of  length  n  followed  by  2  log^p  -  5  real  multi¬ 
plications  for  each  inner  product.  Counting  one  complex  multiplication  as 
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four  real  multiplications  and  considering  the  special  form  of  the  W  col¬ 
umns  gives  a  total  of 

4n(m+I)  +  (2  log^P  -  5)(m+n+JO  (34) 

real  multiplications.  Clearly  the  special  methods  are  substantially  better 
than  the  straightforward  methods  when  p  ^  8  and  m  >  n.  Even  for  p  ■  4, 
the  special  methods  are  superior  because  the  factor  2  log^p  -  5  in  (34) 
must  be  replaced  by  zero.  In  the  derivation  of  both  (33)  and  (34)  it  was 
assumed  that  the  last  row  of  (24)  in  the  dual  problem  was  specially  treated 
to  avoid  multiplications  by  1  and  0.  Also,  (33)  is  valid  for  all  p  ^  3, 
while  (34)  is  valid  only  for  p  =  8,16,32,...  . 

We  earlier  postponed  the  derivation  of  simplex  multipliers 
X(0)  £  R2n+1  corresponding  to  the  initial  basic  feasible  solution  (25). 
They  are  now  derived  in  the  complex  format  (26).  Multiplying  the  initial 
basis  inverse  on  the  left  by  the  row  vector  containing  the  cost  coeffi¬ 
cients  of  the  initial  basic  variables  gives  the  row  vector  .  The 

initial  basis  inverse  is  the  identity  matrix,  the  cost  coefficients  of 
the  basic  W  variables  are  given  in  Table  1,  and  the  cost  coefficient  of 
the  slack  variable  Q  is  0.  Since  8^*0  and  “  ^/2,  we  have 

X(0)  »  (d+(a)R,  d+(-ia)R,  0)  ■  (d+aR,  d+a1,  0)  c  R2n+*  . 

The  definition  (26)  thus  gives 

z(0)  *  a  +  d  eln/4  n  £  C"  ,  £<0)  A  0  .  (35) 
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From  the  proof  of  Theorem  3  it  can  be  seen  that  e  »  0  for  as  long  as  the 
slack  variable  remains  in  the  basis  and  is  positive. 


The  only  sparse  matrices  in  the  dual  problem  are  the  matrices  S,  W, 
and  T.  Any  basic  feasible  solution  of  the  dual  consists  of  2n+l  vari¬ 
ables,  all  of  which  must  be  non-negative.  Since  the  remainder  of  the 
(m+n+A)p  +  1  dual  variables  must  be  zero,  the  matrices  S,  W,  and  T  must  be 
sparse.  Furthermore,  no  row  of  S,  W,  and  T  need  contain  more  than  two 
positive  entries  as  the  next  theorem  shows. 

Theorem  4.  No  basic  feasible  solution  of  the  dual  problem  (23)  -  (24)  can 
have  more  than  two  basic  variables  in  any  one  row  of  W  or  T.  If  a  basic 
feasible  solution  of  the  dual  problem  has  corresponding  simplex  multipliers 
(26)  with  e  >  0,  then  it  cannot  have  more  than  two  basic  variables  in  any 
one  row  of  S. 

Proof .  We  prove  the  first  statement  for  the  matrix  T,  since  the  proof 
for  the  W  matrix  is  just  a  special  case.  Consider  the  j-th  row  of  T. 
Suppose  a  basic  feasible  solution  exists  which  has  the  three  basic  vari¬ 
ables  T.  ,  T.  ,  and  T,  with  r,  s,  and  t  being  distinct.  Then  the 
jr’  js’  jt 

corresponding  simplex  multipliers  z  e  Cn  and  e  e  R  must  result  in  zero 
reduced  costs  for  all  three  variables.  A  result  analogous  to  (27)  was  used 
to  prove  (30);  using  that  analogous  result  here  too  gives 


0 


-i9 

(,Bj '  8i>  •  11 


R 


q  -  r,s,t  . 


(36) 
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Thus  the  single  complex  number  zB^  -  must  have  the  same  projection, 
namely  c^  ,  in  each  of  three  distinct  directions.  This  is  impossible  un¬ 


less  zB^  -  gj  »  0^  “0,  in  contradiction  to  our  initial  assumption  that 
Cj  >  0.  This  establishes  the  first  statement  of  the  theorem.  The  second 
theorem  statement  is  proved  in  the  same  way,  by  using  (27)  itself.  This 
concludes  the  proof. 

The  following  theorem  relates  knowledge  of  an  optimal  basis  of  the 
dual  to  "observable"  quantities  in- the  primal  problem.  The  results  of  the 
theorem  depend  on  the  names ,  but  not  the  actual  numerical  values ,  of  the 
optimal  basic  variables.  In  addition  it  seems  to  indicate  that  the  upper 
bound  CIA)  in  Theorem  2  will  often  be  attained  in  practice. 


Theorem  5.  Let  e'  c  R  and  2'  £  c”  denote  the  simplex  multiplier’s  in 
the  form  (26)  of  a  given  optimal  basis  for  the  dual  problem  (23)  -  (24), 
and  suppose  that  e '  >  0.  If  the  j-th  row  of  one  of  the  matrices  S,  W, 
or  T  contains  two  optimal  basic  variables  in  columns  r  and  t  with 
p  _>  r  >  t  >_  1,  then  either  r-t-1  or  r  -  t  -  p-1.  If  r  -  t  ■  1, 
then 

z'Aj  -  fj  -  e'  sec(u/p)  exp[i(2t  -  l)x/p)  (37) 

z'Bj  -  gj  -  Cj  sec(x/p)  exp[i(2t  -  l)x/p]  (38) 

zj  -  a^  -  dj  sec(x/p)  exp[i(2t  -  l)x/p),  (39) 

according  to  whether  the  j-th  row  is  a  row  of  S,  T,  or  W,  respectively. 

Replacing  t  with  p  in  (37)  -  (39)  gives  the  equations  corresponding  to 
the  alternative  case  r  -  t  -  p  -  1. 
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Proof .  We  treat  only  the  S  matrix  case  since  the  other  two  cases  are 


very  similar.  Thus,  the  two  basic  variables  involved  are 

and  we  assume  that  p  r  >  t  _>  1.  The  reduced  costs  C;Jr 
must  be  0,  so  from  (27) 


S .  and 
jr 

and  c|C 

w 


-i0 


c  ' 


(Z  '  Aj  “  V6 


•~i0 

(-•W 


(40) 


Any  complex  number  having  identical  projections  in  two  directions  is 
uniquely  defined  in  both  magnitude  and  phase.  If  0^  differs  from  0^ 
by  x  radians,  the  system  (40)  implies  that  e’  *  0,  contrary  to  assump¬ 
tion.  Thus  9^  and  9^  differ  by  more  or  less  than  it,  and  (40)  implies 
that 

jz'A^  -  fjj  -  e’  sec(4</2)  >  0  , 


where  <l>  -  min  (0^  -  9^  ,  2x  -  9^  +  0  }.  By  Theorem  2,  4<  must  equal  x/p, 

so  that  either  9  -  0  »  x/p  or  9  -  0  *  x(2p  -  1 ) / p .  From  (11)  we 
rt  rt 

have  either  r  -  t  *  1  or  r-t-p-1.  For  r  -  t  -  1,  solving  the 
system  (40)  for  the  phase  of  z'A^-  f^  gives  (37).  The  case  r  -  t  ■  p  -  1 
is  handled  in  the  same  way.  This  completes  the  proof. 


Theorem  5  is  useful  in  practice  as  well.  Computed  optimal  basic 
solutions  can  be  inspected  very  easily  to  verify  that  the  optimal  basic 


variables  occurring  in  the  same  row  are  in  fact  "paired"  in  the  manner  des¬ 
cribed.  If  they  are  not,  then  one  must  conclude  that  premature  termination 
occurred  in  the  simplex  algorithm,  or  else  that  numerical  round-off  errors 
have  adversely  affected  the  computed  solution. 

Theorem  6.  Let  e'  and  z'  be  as  in  Theorem  5.  If  the  j-th  row  of  one 
of  the  matrices  S,  W,  or  T  contains  an  optimal  basis  variable  in  column 
r,  1  <_  r  _<  p,  then 

e'  *  b’Aj  -  fjUe’ sec 

9r  *  */p  £  arS  <2 1  a  -  f j)  £  0r  +  it/p,  (41) 

or 

-  li,BJ  - 

9r  *  <  ar8  (z'Bj  -  gj)  <  9r  +  n/ p,  (42) 

or 

dj  <_  |zj  -  aj  j  <  dj  sec  it/p 

9r  "  ny,P  1  ar8  (zj  "  aj)  <.  9r  +  */p.  (A3) 

according  to  whether  the  j-th  row  is  a  row  of  S,  W,  or  T,  respectively. 

Proof .  The  proof  is  closely  related  to  the  method  of  proof  of  Theorem  5, 
and  it  is  not  presented.  We  remark  that  if  Theorem  6  were  proved  first, 
then  Theorem  5  would  be  an  Immediate  consequence  of  it. 

This  section  is  concluded  with  an  alternative  statement  of  the  dual 
problem  (23)  -  (24)  using  complex  arithmetic  notation.  We  do  not  use  it 
elsewhere  in  this  report,  but  this  form  is  interesting  in  its  own  right  and 
also  provides  a  more  concise  statement  of  the  dual  which  may  be  useful. 
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Dual  Problem:  complex  format 


min 

S,T 

W,Q 


-10  R 

[(fS+gT+aW)e  ]  + 


f  (CT  +  DM  ) 
j  =  l  J  J 


subject  to:  S  ^  0,  T  ^  0,  W  ^  0,  Q  ^  0,  and 

(AS  +  BT  +  W)e“i9  -  0  e  C° 
m 

Q+  I 

j-l  k-1 


1  £  R. 


We  have  used  e  to  denote  a  complex  column  vector  of  length  n  whose 

“i9k 

k-th  component  is  e  e  C;  all  other  notation  is  unchanged  from  (23)  - 
(24).  Verifying  that  this  complex  format  is  correct  is  straightforward. 

III.  Solution  of  the  linearized  problem  for  large  p. 

The  relative  accuracy  between  the  absolute  value  and  the  linearized 
absolute  value  is  given  by  (12).  Although  in  some  applications  a  choice  of 
p  *  16  may  provide  satisfactory  accuracy,  in  other  applications 
p  -  16  may  be  much  too  small.  As  stated  in  the  Introduction,  p  *  1024 
gives  5  significant  digits  of  relative  accuracy,  so  applications  re¬ 
quiring  greater  accuracy  will  therefore  need  p  >  1024.  There  is,  however, 
a  practical  limit  to  how  large  p  may  be  taken  in  many  problems.  For  ex¬ 
ample,  a  problem  will  become  numerically  unstable  for  sufficiently  large 
values  of  p  if  its  optimal  solution  has,  for  every  p,  two  basic  dual 
variables  in  at  least  one  row  of  S  W,  or  T.  From  Theorem  5,  these  prob¬ 
lems  have  two  "consecutive''  projections  in  an  optimal  basis  and,  since  the 
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inequalities  defining  these  projections  become  progressively  less  distin¬ 
guishable  numerically  as  p  increases,  the  basis  matrix  must  become  pro¬ 
gressively  more  ill-conditioned.  Only  those  problems  which  never,  for  any 
p,  have  more  than  one  optimal  basic  variable  in  any  one  row  of  S,  W,  or  T 
will  escape  numerical  ill-conditioning  from  this  cause;  however,  problems 
of  this  kind  seem  to  be  uncommon. 

It  is  useful,  nonetheless,  to  be  able  to  solve  the  linearized  problem 
for  as  large  a  value  of  p  as  is  numerically  practical.  One  reason  is 
that  a  solution  of  the  linearized  problem  furnishes  a  starting  point  for 
other  methods  which  potentially  provide  greater  accuracy.  For  instance, 
the  problem  (1)  -  (4)  can  be  rewritten  as  a  semi -inf inite  program,  or  SIP, 

and  an  interesting  algorithm  [10],  [11)  for  solving  a  class  of  general 

• 

SIP's  can  be  utilized  to  solve  it.  This  method  sets  up  an  appropriate 
nonlinear  system  of  algebraic  equations  which  are  then  solved  using  the 
Newton-Raphson  method  (or  other  workable  iterative  method).  A  feasible 
solution  of  this  nonlinear  system  gives  a  solution  of  the  SIP.  The  start¬ 
ing  point  of  the  Newton-Raphson  iteration  can  be  taken  to  be  the  solution 
of  the  linearized  problem  (7)  -  (10).  Conceivably  a  very  good  starting 
point  may  be  needed  to  ensure  that  the  Newton-Raphson  iteration  converges 
to  a  feasible  point.  In  this  event,  the  linearized  problem  can  be  made 
more  accurate  merely  be  Increasing  p.  The  subproblem  (1)  -  (2)  is  treated 
as  an  SIP  in  [5],  alghough  no  mention  is  made  there  of  special  structure  in 
the  linearized  problem  or  of  results  such  as  those  in  any  of  the  theorems 
above. 
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The  method  we  suggest  for  solving  the  linearized  problem  for  large 
values  of  p  begins  by  solving  the  smallest  dual  problem.  That  is,  the 
dual  problem  (23)  -  (24)  with  p  *  4  is  solved  by  the  simplex  method 
starting  at  the  basic  feasible  solution  (25).  Next,  the  p  ■  8  dual 
problem  is  solved  using  the  optimal  basis  for  the  p  *  4  dual  to  start  the 
simplex  algorithm.  The  p  ■  16  dual  is  then  solved  starting  at  the  optim¬ 
al  basis  for  the  p  »  8  dual.  Continuing  in  this  fashion,  large  values  of 
p  are  quickly  reached.  The  algorithm  is  always  well-defined  because  basic 
feasible  solutions  of  the  dual  for  a  given  p  are  also  basic  feasible 
solutions  of  the  dual  for  all  larger  values  of  p.  This  special  feature 
of  the  dual  holds  because  the  ©  sets  are  nested  for  the  allowed  values 
p  -  4 ,  8 ,  16 ,  32 ... .  . 

By  doubling  p  at  each  stage  beginning  with  p  »  4,  this  algorithm 
permits  the  simplex  iterations  to  avoid  bases  associated  with  numerical  in¬ 
stability  from  the  linearization  process  until  p  becomes  large.  In  other 
words,  potentially  troublesome  bases  are  not  encountered  until  a  signifi¬ 
cant  number  of  simplex  iterations  have  already  transpired.  Therefore 
greater  numerical  accuracy  might  be  anticipated.  The  algorithm,  of  course, 
cannot  avoid  difficulties  caused  by  ill-conditioning  in  the  complex  equa¬ 
tions  themselves. 

This  method  is  practical  for  large  values  of  p  only  because  the 
methods  presented  in  section  II  require  computer  storage  that  is  almost  in¬ 
dependent  of  p.  As  pointed  out  there,  the  dual  problem  requires  only 

2n(m+i)  +  2p  storage  locations.  With  an  explicit  inverse  formulation  of 

2 

the  simplex  method,  an  additional  (2n+l)  words  are  required.  Thus  the 
only  storage  dependent  on  p  are  the  2p  locations  for  storing  the  cosines 


and  sines  of  the  projections.  (Actually,  only  p/8  storage  locations  can 
be  made  to  suffice  for  this  purpose.) 

One  advantage  of  this  algorithm  is  that  the  optimal  basis  for  each 
intermediate  value  of  p  can  be  easily  inspected  using  Theorems  5  and  6  to 
determine  to  what  extent  numerical  round-off  errors  are  present  in  the 
solution.  If  sufficient  error  is  present,  the  algorithm  can  be  terminated 
early,  or  alternatively,  the  basis  can  be  reinverted  before  continuing  to 
the  next  value  of  p.  Either  way,  pointless  simplex  iterations  can  be 
avoided. 

The  only  drawback  the  algorithm  potentially  has  is  that  more  simplex 
iterations  might  be  required  to  reach  the  final  optimal  dual  basis  by 
forcing  it  to  proceed  via  smaller  values  of  p  than  by  solving  the  full 
dual  problem  all  at  once.  This  difficulty  does  not  seem  to  be  significant 
in  practice.  Should  it  ever  become  a  problem,  however,  it  could  be  at 
least  partially  overcome  by  skipping  more  rapidly  through  the  available 
values  of  p.  For  instance,  rather  than  doubling  p  at  each  stage,  one 
could  quadruple  it  instead.  It  is  also  possible  to  begin  the  algorithm 

with  a  larger  initial  value  of  p;  that  is,  p  >  4, 

An  incidental  advantage  of  the  algorithm  is  that  the  primal  solutions 
for  all  the  intermediate  values  of  p  can  be  stored.  This  raises  the  pos¬ 
sibility  of  solving  the  original  problem  (I)  -  (4)  by  extrapo1 ation.  The 
original  problem  corresponds  to  the  case  p  -  ®,  so  it  is  not  unreasonable 
to  think  that  better  solutions  could  be  obtained  by  extrapolating  the 

solutions  corresponding  to  the  finite  cases  p  -  4,  8,  16,  32,...  .  The 

optimal  vectors  z  of  the  linearized  primal  problems  converge  linearly 
with  increasing  p,  while  the  optimal  values  e  converge  quadratically 
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with  increasing  p.  Therefore  Richardson  extrapolation  (see,  for  instance 
[12]  or  [13]) might  be  suitable  for  extrapolating  primal  solutions.  Extrap¬ 


olating  the  p  =  ®  dual  solution  [14]  from  the  finite  p  dual  solutions 
requires  special  care.  The  optimal  finite  dual  basis  names  must  be  normal¬ 
ized  before  extrapolation,  i.e.,  the  column  numbers  of  the  basic  variables 
must  be  multiplied  by  2x/p  to  map  them  to  the  interval  (0,2x),  and  the 
possible  discontinuity  at  2u  must  be  taken  into  account  during  extrapola¬ 
tion.  Also,  those  optimal  finite  p  dual  basic  variables  which  occur  in 
"pairs",  as  described  in  Theorem  5  above,  must  coalesce  into  one  optimal 
p  ■  ®  dual  basic  variable. 

IV.  Details  of  the  revised  simplex  algorithm 

It  is  not  desirable  to  use  a  computer  code  which  treats  the  complex 

0 

matrices  and  vectors  of  the  primal  problem  by  separating  them  into  their 
real  and  imaginary  parts.  Besides  the  inconvenience,  such  a  code  would  be 
very  inefficient  in  practice  because  the  reduced  cost  coefficient  calcula¬ 
tions  (29)  -  (30)  would  cause  thrashing  on  virtual  memory  systems.  (See 
[15]  for  an  example  involving  Hermitian  eigenproblem  solution.)  Given, 
then,  that  the  solution  vector  z  of  the  primal  problem  is  best  stored  as 
a  complex  vector,  it  becomes  clear  that  the  simplex  multipliers  (26)  should 
be  re-ordered  to  reflect  the  storage  of  z.  Consequently  the  rows  (24)  of 
the  dual  problem  should  also  be  re-ordered.  The  computer  code  therefore 

visualizes  the  dual  problem  rows  in  the  following  order: 

* 

{ 1 ,  n+1,  2,  n+2,  ...,  n-1,  2n-l,  n,  2n,  2n+l}  (44) 
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where  these  numbers  denote  the  row  numbers  in  the  original  system  (24). 

The  re-ordered  system  turns  out  to  be  much  easier  to  work  with  in  FORTRAN 
code  than  the  original  system.  With  the  rows  of  the  dual  problem  ordered 
as  in  (44),  the  reduced  cost  calculations  can  be  coded  in  FORTRAN  just  as 
they  are  written  in  (29)  -  (31),  provided  the  initial  data  of  the  problem 
are  typed  COMPLEX. 

The  most  negative  reduced  cost  (32)  determines  the  entering  basic 
variable  in  the  simplex  algorithm.  If  no  reduced  cost  is  negative,  the 
simplex  algorithm  terminates  because  the  current  basis  is  optimal.  If  ties 
develop  for  the  most  negative  reduced  cost,  they  are  broken  by  choosing 
the  variable  with  the  least  lexicographical  index.  Since  every  dual  vari¬ 
able  has  three  names,  the  "index"  of  a  dual  variable  is  a  triplet  of  posi- 

tive  integers: 

i/j/k 

where  i  *  1,  2,  or  3  according  to  whether  it  is  an  S,  W,  or  T  variable 

j  »  constraint  number;  from  (8)  -  (10) 

k  »  projection  number  of  the  angle  in  the  set  0;  1  £  k  £  p. 

Note  that  the  middle  index  j  has  different  ranges  of  possible  values  de¬ 
pending  on  the  value  of  i.  These  triplets  are  ordered  lexicographically, 
so  the  least  index  is  well-defined. 

There  is  one  exception  to  the  least  index  rule  in  case  of  ties  for 
the  entering  variable.  The  exception  arises  because  the  highly  degenerate 
initial  starting  point  (25)  can  cause  cycling  in  the  simplex  algorithm.  As 
long  as  the  slack  variable  Q  remains  in  the  basis,  the  only  entering 
variables  permitted  are  S  variables  with  negative  reduced  costs.  If  S 
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variables  with  negative  reduced  costs  do  not  exist,  then  the  entering  vari¬ 
able  is  permitted  to  be  a  W  or  a  T  variable  and  the  ties  are  resolved 
by  the  least  index  rule  as  described  above.  Thus,  S  variables  are  given 

priority  for  entering  the  basis  only  for  as  long  as  the  slack  Q  is  in  the 
basis.  Once  Q  is  removed  from  the  basis  it  is  never  allowed  to  re-enter, 
and  the  exception  to  the  least  index  tie  breaking  rule  becomes  moot. 

The  outgoing  basic  variable  is  determined  by  the  usual  ratio  test, 
with  ties  resolved  by  taking  the  variable  having  the  maximum  magnitude 
pivot  with  the  smallest  index.  That  is,  the  variable  having  the  least 
ratio  in  the  ratio  test  must  leave  the  basis.  If  the  least  ratio  is  at- 
tained  by  more  than  one  variable,  the  variable  having  the  largest  magnitude 
pivot  is  chosen.  If  more  than  one  of  these  variables  have  the  same  magni¬ 
tude  pivot,  then  the  variable  with  least  index  is  selected.  . 

Because  of  degeneracy  and  cycling,  there  is  one  exception  to  this  tie 

breaking  rule  for  the  exiting  variable.  So  long  as  the  slack  Q  remains 
in  the  basis,  only  W  variables  are  permitted  to  exit.  This  rule  makes 

sense  only  when  a  W  variable  is  involved  in  the  tie;  if  no  such  W  vari¬ 
able  exists,  the  exception  is  not  invoked.  If  more  than  one  W  variable 
is  involved  in  the  tie,  then  the  one  having  the  largest  magnitude  pivot 
with  the  least  index  is  selected  to  exit.  Just  as  for  the  entering  vari¬ 
able,  this  exception  becomes  moot  once  the  slack  Q  leaves  the  basis. 

With  these  modifications  to  the  usual  tie  breaking  rules  for  entering 
and  exiting  variables,  no  cycling  in  the  simplex  algorithm  has  yet  been  ob- 
served.  However,  if  these  modifications  are  deleted,  cycling  may  very  well 
occur.  Example  3  of  section  V  below  cycled  (with  a  cycle  of  length  19) 
without  these  modifications.  Whether  or  not  cycling  in  this  particular 


example  still  happens  when  exact  arithmetic  is  used  is  not  known  to  the 
author.  It  is  possible  that  the  observed  cycling  is  merely  an  artifact  of 
finite  precision  arithmetic. 

In  practical  implementations  of  the  simplex  algorithm,  a  zero 
tolerance  is  necessary  when  testing  for  the  most  negative  reduced  cost  and 
for  possible  divisors  in  the  ratio  test.  This  number  is  somewhat  arbi¬ 
trary,  but  it  is  important  that  it  not  be  too  small  and  that  it  somehow  be 


dependent  on  the  scale  of  the  data  of  the  problem.  The  number  used  here  is 
the  product  of  the  unit  round-off  error  of  the  host  computer  and  the  sum  of 
the  absolute  values  of  the  incoming  column  (i.e.,  its  norm).  This 
number  is  rather  conservative  but  seems  satisfactory  in  the  examples  run  to 
date.  It  is  used  for  both  the  reduced  cost  and  the  pivot  tolerance  tests. 

Besides  the  usual  expected  termination  criteria  in  the  simplex  al¬ 


gorithm,  the  pricing  method  implicit  in  (29)  -  (31)  yields  an  interesting 


and  perhaps  novel  way  to  terminate  the  algorithm.  The  pricing  method  com¬ 
putes  the  most  negative  reduced  cost  by  examination  of  all  the  reduced 
costs,  not  merely  the  reduced  costs  of  the  nonbasic  variables.  Hence  it 


can  happen  that  the  entering  and  the  exiting  variable  are  identical  because 
of  numerical  round-off  errors.  In  practice  this  event  seems  to  signal  that 
no  further  improvement  in  the  solution  is  numerically  possible.  Solutions 
returned  by  terminating  the  algorithm  whenever  this  "self-cycling"  occurs 
appear  to  be  very  satisfactory.  On  the  other  hand,  failure  to  terminate  in 


this  situation  leads  to  wasted  simplex  iterations.  Eventually  enough  ac¬ 
cumulated  numerical  round-off  error  develops  to  end  self-cycling,  but  by 
then  the  solution  is  usually  very  poor. 
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A  FORTRAN  code  was  developed  to  implement  and  test  the  methods  de¬ 


scribed  for  solving  the  dual  problem.  It  holds  an  explicit  basis  inverse 
and  performs  the  usual  pivoting  to  update  the  inverse  in  each  simplex 
iteration.  This  procedure  is  known  to  be  numerically  unstable,  but  easily 
programmed.  To  forestall  numerical  difficulties  the  inverse  is  held  in 
double  precision,  even  though  a  double  precision  inverse  is  not  in  general 
a  satisfactory  substitute  for  a  numerically  stable  technique.  Updating  the 
LU  or  QR  factorizations  of  the  basis  are  preferable.  Nonetheless  the  ex¬ 
plicit  inverse  code  gave  good  performance  on  the  small  problems  tested  in 
this  report.  Reinversion  of  the  basis  is  also  a  desirable  feature  of  a 
general  purpose  code,  but  it  was  not  made  part  of  the  initial  FORTRAN  code. 


V.  Examples 

'  0 

Example  1  is  taken  directly  from  [4,  p.  249].  It  can  be  put  in  the 

form  (1)  -  (3)  by  letting  n  *  2,  m  »  5,  X  *  2,  and  defining  the  matrices 


f 


-1  +  i 
-1  +  i 


.5i 

0 


(45) 


-1  +  i 


In  this  example  all  the  data  happen  to  be  real  valued,  except  the  vector 
f.  Since  constraints  of  type  (4)  are  missing  from  this  example  and  since 
their  linearizations  provide  the  initial  dual  basis,  we  must  insert  them 
here.  We  take 


a 


0 


(46) 


so  that  they  will  be  inactive  at  the  solution.  The  exact  solution  of  this 
problem  is  z^  ■  (-1  +  i)/2,  z^  “  0,  and  e  *  HU 

Table  2  gives  the  solutions  of  the  linearized  primal  problem  for 
p  »  4,  8,  16,..., 2048.  It  is  apparent  that  the  optimal  value  of  e  for 

p  =  8  is  also  the  optimal  value  of  z  for  all  p  16.  For  p  _>  8,  then, 

the  accuracy  of  the  primal  solutions  depends  solely  upon  the  linearization 

errors  since  the  optimal  e  does  not  change.  Furthermore,  the  optimal 

value  of  the  original  nonlinear  problem  must  surely  equal  the  optimal  z 
for  p  »  8. 

Table  3  gives  the  optimal  basic  solutions  of  the  dual  problems  for 

p  *  4,  8,..., 2048.  The  active  constraints  are  the  same  for  all  p  >_  8, 

except  for  their  9  names.  Hence  the  active  constraints  at  the  optimum 
of  the  original  nonlinear  problem  (1)  -  (4)  have  surely  been  identified. 
The  fourth  and  fifth  basic  variables  are  "paired"  in  an  obvious  way,  this 
behavior  is  explained  by  Theorem  5. 
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All  the  optimal  solutions  presented  in  Table  3  are  degenerate,  or 
very  nearly  so.  It  is  particularly  interesting  that  the  "degenerate  parts" 
of  the  optimal  solutions  almost  double  in  size  as  p  is  doubled, 
especially  when  p  >_  64.  Assuming  this  trend  continues,  the  optimal 
solution  will  eventually  look  as  if  it  is  not  degenerate.  We  suspect,  but 
have  not  attempted  to  prove,  that  this  trend  is  purely  an  artifact  of  the 
numerical  ill-conditioning  inherent  in  the  linearization  process,  as  was 
discussed  in  section  III  above. 

Assuming  linear  convergence  for  increasing  p,  one  step  of  the 
Richardson  extrapolation  process  is  easily  applied  to  the  last  two  z 
vectors  in  Table  2.  Multiplying  the  p  »  2048  vector  by  two  and  sub¬ 
tracting  the  p  »  1024  vector  gives 


R 

Z1 

-.500000 

I 

Z1 

+.500C00 

R 

Z2 

+.001370  *  10_L1 

I 

Z2 

i  4 

1 

o 

.  4 

X 

00 

o 

o 

o 

• 

1 

__  __ 

L  J 

Evidently  one  step  of  this  extrapolation  procedure  has  nearly  doubled  the 
number  of  correct  significant  digits  in  the  primal  solution. 
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It  is  possible  to  make  this  problem  infeasible  by  adjoining  only  one 


constraint  of  type  (3).  Instead  of  (45)  we  take 


2  2 

1 

.  g  * 

0 

,  c  = 

/2 

2  -4 

1 

0 

7-4i 

/2 

29/4 

The  linearized  problem  is  feasible  only  for  p  -  4  and  8;  for  p  >_  16,  it  is 
infeasible.  This  illustrates  the  repark  made  in  the  Introduction  that  some 
linearized  problems  may  have  feasible  solutions  even  when  the  original 
problem  (1)  -  (4)  is  actually  infeasible. 

Example  2  is  the  same  as  Example  1,  except  that  constraints  of  type 

(4)  are  tightened  so  that  they  are  active  at  the  solution.  Instead  of 

# 

(46),  we  take 


a  ■ 

1 

o 

.  d  - 

1 

• 

0 

.4 

—  — 

_  _ 

The  exact  solution  of  this  problem  is  e 


✓2  -  .4,  z i  -  (-1  +  i)  /I/5 


and 

1  _  317  (/2  -  1)  -  (431902  -  190320/1) 1/2 

2  3000  -  1200/1 


-  -.208846903 


zR  -  -  'LL  +  [  i  -  (z1  -  LL  )  ]  -  -.093336568 

2  10  8  2  10 


Tables  4  and  5  give,  respectively,  the  solutions  of  ehe  linearized  primal 


and  dual  problems  for  p  »  4,  8,...,  2048.  As  in  Example  1,  the  optimal  e 
for  p  *  8  is  optimal  also  for  all  p  _>  16.  The  Richardson  extrapolation  of 
the  last  two  z  vectors  in  Table  4  can  be  performed  as  described  in  Ex¬ 
ample  1  above.  In  this  case  one  extrapolation  step  gives 
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R 
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i 

zi 
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+.282843 

R 

Z2 

+.0933341 

I 

Z2 

-.208848 

_  _ 

_  _ 

which  is  an  improvement  of  2  to  3  significant  digits  over  the  best  answer 


in  Table  4 


ON 

o 


GO 

o 

CM 


r- 

r**. 

cm 

eo 

00 

CM 


CM 

GO 

CM 


ON 

-O 

O 

CO 

CM 

ON 

O 


o 

sO 

ON 

00 

o 

CM 


• 

• 

• 

• 

• 

1 

1 

1 

*H 

r«. 

o 

m 

in 

CM 

H 

h- 

r-. 

CM 

ON 

CM 

O 

CM 

o 

o 

CO 

ph 

CM 

ON 

CO 

00 

00 

ON 

o 

ph 

CM 

CM 

O 

CM 

o 

• 

• 

• 

• 

• 

1 

i 

1 

pH 

in 

CM 

CO 

CO 

O 

•o 

nO 

CM 

-4- 

ON 

fH 

CM 

uO 

pH 

CM 

CM 

CM 

NO 

<r 

pH 

ON 

<• 

CM 

m 

00 

00 

ON 

o 

>—* 

CM 

CM 

o 

CM 

O 

• 

• 

• 

• 

• 

1 

1 

1 

pH 

00 

CM 

NO 

CO 

<r 

*—♦ 

r>* 

CO 

CM 

t—i 

v£> 

co 

CO 

CM 

in 

m 

no 

ON 

ON 

ON 

CM 

CM 

00 

r*^. 

00 

O 

pH 

CM 

CM 

o 

CM 

o 

• 

• 

• 

• 

• 

1 

1 

f 

H 

CO 

NO 

ON 

pH 

CM 

*3- 

00 

ON 

vO 

00 

pH 

oo 

r** 

CO 

0> 

in 

CM 

pH 

CM 

ON 

in 

^3“ 

o 

<r 

CM 

H 

00 

Cv 

CO 

1—* 

ph 

CM 

CM 

O 

CM 

o 

• 

• 

• 

• 

• 

1 

1 

i 

ph 

ON 

00 

00 

o 

CM 

<r 

co 

NO 

nO 

ph 

<r 

r^« 

ON 

CM 

OO 

CM 

00 

nO 

vO 

00 

pH 

<r 

pH 

ON 

NO 

pH 

h 

CM 

CM 

o 

IN 

o 

• 

• 

• 

• 

• 

i 

1 

i 

H 

«H 

o 

m 

r-* 

O 

•4- 

o 

00 

m 

ON 

pH 

CM 

r- 

ON 

CO 

CO 

CM 

NO 

<o 

O 

*o 

H 

<r 

<r 

pH 

pH 

m 

nO 

pH 

H 

CO 

CM 

o 

CM 

o 

• 

« 

• 

• 

• 

1 

1 

1 

pH 

H 

»H 

00 

ON 

ON 

o 

H 

vO 

r- 

CO 

in 

NO 

CM 

CO 

H 

CO 

H 

nO 

pH 

ON 

r-* 

H 

pH 

CM 

CM 

O 

CM 

o 

• 

• 

• 

• 

• 

;  • 

i 

i 

! 

CM 

CM 

H 

-3- 

CO 

fH 

CM 

CM 

CM 

pH 

00 

m 

o 

NO 

CO 

H* 

H 

<r 

CM 

CM 

*3- 

pH 

CO 

CM 

o 

CM 

o 

• 

• 

• 

• 

• 

i 

1 

o 

o 

CO 

CO 

© 

o 

o 

in 

in 

o 

o 

o 

m 

m 

o 

o 

o 

CO 

CO 

o 

n- 

o 

o 

m 

m 

o 

H 

H 

NO 

• 

• 

• 

1 

i 

CM 

•O 


CO 

CO 


at  h 

N 


♦H  0$  CM 

N  N 


W  CM 
N 


I 

•H  CO  (B 

fl  w  c 
u  V  o 
O  «H 
u)  H  m  w 


o. 

e 

cG 

X 

CxJ 


s 

0J 

.O 

o 

a 


CO 

B 

«H 

Ui 

a 

■o 

o> 

M 

«H 

U 

CO 

a 

c 


QJ 

-C 


CO 

c 

o 


o 

CO 


■o 

0) 

pH 

Si 

CO 

H 


43 


Example  3  is  taken  from  [6J  and  is  an  unconstrained  Chebyshev  func¬ 


tion  approximation  problem;  that  is,  constraints  (3)  -  (4)  are  absent.  The 
101  columns  of  the  matrix  A  e  c^*^l  are 


1 

exp[i(j-l)*/4001 

exp[i2(j-l)?t/400] 


»  j  “  1,2, ...,101 


while  the  components  of  f  e  ClOl  are 


fj  -  exp  [i3( j-l)x/400]  ,  j  -  1,2, ...,101. 

i3x 

In  other  words,  the  complex  valued  function  e  is  approximated  by  com¬ 
plex  linear  combinations  of  the  three  functions  { 1 ,e^x,e^  x}  over  101 
equispaced  points  on  the  x-interval  [0,n/4].  As  in  Example  1  we  insert 
bounds  of  type  (4)  to  provide  an  initial  dual  basis;  specifically,  we  .take 


—  — 

—  — 

0 

.  d  - 

10 

0 

10 

0 

10 

—  — 

_ 

These  constraints  are  not  active  at  the  optimal  solution.  It  is  not 

necessary  to  insert  artificial  inequalities  of  type  (3). 

It  can  be  verified  that  the  explicit  solution  of  Example  3  can  be 
written 

z1  -  Tj  exp  (i  3x/8) 
z2  -  r2  exp  (i  5it/4) 
z3  -  rj  exp  (i  n/8) 
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I 


where 


r  =■  a  =  .96157056080646 

r2  =  b-2(b-a2)/(l-a2)  =  2.8122548927058 

r3  -  a(l-2b+a2)/( 1-a2)  =  2.8477590650226 

a  =*  Xcos(u/ 16)+(  l-X)cos(it/8) 
b  =  Xcos(u/8)+( l-X)cos(x/4) 
c  ■  Xcos(3u/ 16)+( l-X)cos(3x/ 8) 

X  =  (sin(ir/8))/(sin(x/16)+sin(x/8)) 

e  =  ( 1-cr L+  br2  -  ar3)1/2  -  .014706309694449 

Tables  6  and  7  give,  respectively,  the  solutions  of  the  linearized 
primal  and  dual  problems  for  p  -  4,  8,  ...,  1024.  The  optimal  value  of  e 
for  p  -  64  seems  to  be  the  optimal  e  value  for  all  p  >  128.  The  obvious 
"pairing"  of  the  basic  variables  in  Table  3  is  explained  by  Theorem  5. 

Note  also  that  the  basis  values  are  not  altered,  only  permuted,  for 
p  >  32. 
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Richardson  extrapolation  can  be  applied  to  the  last  two  z  vectors 
in  Table  6  just  as  in  the  two  previous  examples.  In  this  case  one  extrapo¬ 
lation  step  gives 


— 

—  — 

R 

21 

.367978 

I 

21 

.888376 

R 

Z2 

m 

-1.988564 

I 

Z2 

-1.988564 

R 

z3 

2.630986 

I 

Z3 

1.089791 

__  _ 

_  _ 

so  that 

jzj  -  .96157148 

| z 2 1  -  2.81225418 
| z 3 1  -  2.84775908  . 

It  is  easy  to  verify  that  before  extrapolation  the  case  p  *  1024  gives 


JzJ  -  .96236 

| z 2 1  -  2.81376 
Jz3j  -  2.84855  . 


Extrapolation  in  this  case  more  than  doubled  the  number  of  correct  signifi¬ 
cant  digits  in  the  solution  vector  z. 

Numerical  computations  for  the  above  examples  were  performed  on  a 
DEC  10  at  Stanford  University. 
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Table  7.  Optimal  solutions  of  the  dual  problem  of  Example 


Another  unconstrained  Chebyshev  function  approximation  problem 
given  in  [6]  is  moderately  large  and  100%  dense.  In  this  example  the 

44x501 

501  columns  of  the  matrix  A  e  C  are  given  by 


Aj  - 

exp  (i  k1  Xj) 
exp  (i  lt2  x  ) 

-  exp  (i  k45  Xj ) 

1 

1 

• 

• 

• 

exp  (i  kA4  x.) 

—  _ 

• 

• 

_1_ 

,  j 


1,2,. ..,501 


where  1  =  k  <  k„  <  . . .  <  k,  ,  <  k,  _  a  49  are  the  distinct  integers  be- 
1  2  44  45 

tween  1  and  49,  excluding  the  integers  {7,17,21,29},  and  where 
X1  “  U0  +  u0)/250*  i  *  1,2,... ,501  with  uQ  -  .0538117. 

501 

The  components  of  f  e  C  are  f^  *  exp  (i  k^,.  x ^  )  ,  j  **  1,...,501. 

This  example  lacks  constraints  of  type  (3)  -  (4).  The  linearized  problem 
for  p  ••  16  was  solved  on  a  DEC  VAX  11/780  in  1350  simplex  iterations. 

Total  CPU  time  was  25  minutes  and  .7  million  page  faults  were  Incurred. 

Only  80,000  words  of  storage  were  needed.  In  contrast,  the  algorithm  pro¬ 
posed  in  [6]  (which  utilizes  the  algorithm  [16]  as  a  subroutine)  solved  the 
same  problem  on  the  same  VAX  in  1270  simplex  iterations,  requiring  179 
minutes  of  CPU  time  and  incurring  11  million  page  faults.  Over  360,000 
words  of  storage  were  needed.  Both  solutions  were  essentially  identical, 
as  expected.  The  difference  in  the  number  of  simplex  iterations  is  ex¬ 
plained  as  follows.  The  algorithm  [6]  solves  the  full  problem  for  p  -  16, 

while  the  algorithm  developed  here  solves  the  p  *  4  problem  and  the  p  ■  8 
problem  before  solving  the  p  ■  16  problem.  This  indirect  route  to  the  full 

problem  solution  is  less  efficient  in  this  example  than  solving  the  full 
problem  immediately. 


VI.  Concluding  remarks 

Extrapolation  was  suggested  as  one  method  for  extracting  more 
accuracy  from  the  linearized  problem  solutions.  An  alternative  procedure 
can  provide  as  much  accuracy  in  the  solution  as  desired,  although  with  more 
work  than  is  required  by  extrapolation.  The  optimal  solutions  of  the 
linearized  dual  problems  for  sufficiently  large  p  must  identify  the  con¬ 
straints  active  at  the  optimal  solution  of  the  original  (i.e.,  p  «  ®) 
problem  (1)  -  (4).  Deleting  the  inactive  constraints  from  the  p  »  ® 
problem  yields  an  equality  constrained  optimiztion  problem  which  the  method 
of  Lagrange  multipliers  seems  well  suited  to  solve.  Lagrange's  method 
gives  rise  to  a  nonlinear  system  of  algebraic  equations  in  the  optimum 
value  e,  the  solution  vector  z,  and  the  multipliers  A.  Iterative  methods 
for  the  solution  of  this  system  can  be  started  from  a  excellent  initial 
point  (e,  z,  A)  provided  by  a  linearized  primal  and  dual  problem  solution. 
Safeguarded  Newton-Raphson  iteration  should  be  a  highly  effective  iterative 
method  for  solving  this  system,  especially  if  advantage  is  taken  of  the 
system's  special  form  (i.e.,  for  A  given,  the  vector  z  can  be  found  by 
solving  a  system  of  linear  equations).  A  possible  limitation  of  this 
approach  is  that  very  large  values  of  p  might  be  necessary  in  order  to 
identify  the  right  p  ■  ®  active  set.  The  examples  of  the  previous  section, 
however,  indicate  that  the  optimal  active  set  is  found  for  relatively  small 
values  of  p.  Specifically,  in  Examples  1,  2,  and  3,  the  correct  p  «  ® 
active  sets  (determined  from  the  optimal  dual  basis  names  in  Tables  3,  5, 
and  7)  first  appear  when  p  is  8,  8,  and  32,  respectively. 

Certain  kinds  of  domain  and  range  constraints  can  be  adjoined  to  the 
constraints  (8)-(10)  with  only  minor  extension  of  the  algorithm  proposed 
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for  solving  (7)  -  (10).  Let  the  matrix  H  «  Cn  q  ,  and  the  row  vectors 
e  £  o'*  ,  £  R**  ,  and  h  c  R^  be  given.  Then  the  constraints 


Re  ((zH..  -  e^  )  exp  (-i^ ))  <  h^  ,  j  -  l,...,q 


R  I 

are  linear  in  z  and  z  ,  and  so  can  be  added  to  the  problem  (7)-(10). 

It  is  easy  to  see  that  the  constraints  (9)  and  (10)  are  instances  of  (47). 
However,  (47)  can  impose  constraints  not  possible  with  (9)  and  (10).  For 
instance,  if  q  ■  1,  the  constraint  that  the  complex  number  zH^-  e^  must 
be  in  the  right  half  of  the  complex  plane  is  equivalent  to 

Re  ( (zH^  -  e^)  exp(-ix))  <  0  .  (48) 

Furthermore,  if  all  the  columns  of  the  matrix  H  are  identical  to  its 
first  column  ,  then  the  number  zH^  -  e^  can  be  confined  to  any  closed 
convex  polygonal  region  (bounded  or  unbounded)  in  the  complex  plane  by 
appropriate  choices  of  h,  and  q  in  (47).  More  generally,  the  con¬ 
straints  (47)  can  be  thought  of  as  confining  each  of  the  complex  numbers 
zHj  -  e^  to  different  closed  convex  polygonal  regions  in  the  plane. 

When  complex  function  approximation  on  a  discretized  arc  or  domain 
boundary  in  the  complex  plane  gives  rise  to  the  problem  (l)-(4),  then  an 
implicit  natural  ordering  of  the  columns  of  the  matrix  A  exists.  This 
ordering  is  inherited  from  the  ordering  of  the  discrete  points  along  the 
arc.  Therefore  it  is  possible  to  devise  clever  strategies  of  both  multiple 
and  partial  pricing  in  the  simplex  algorithm  for  solving  the  linearized 
problem  in  order  to  significantly  reduce  overall  computation  time.  Effect¬ 
ive  partial  pricing  schemes  would  require  far  fewer  evaluations  of  the 
vector-matrix  products  zA^  in  (28)  and  yet  not  incur  significant  in¬ 
creases  in  the  total  number  of  simplex  iteratios  needed  to  reach  opti- 
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mality.  Effective  multiple  pricing  schemes  would  decrease  the  number  of 
simplex  iterations  needed  to  reach  optimality  by  increasing  the  change  in 
e  in  each  simplex  iteration  step.  Both  multiple  and  partial  pricing  can 
be  implemented  simultaneously  and  might  significantly  improve  computation 
time,  especially  when  m  and  n  are  large. 

It  has  been  assumed  throughout  this  report  that  the  unknown  vector  z 
must  lie  in  Cn  .  In  some  applications  it  is  desirable  to  restrict  z  to 
Rn  ,  while  still  retaining  complex  matrices  A  and  B  in  original  problem 
(1)  -  (A).  It  is  easy  to  see  that  setting  z1  -  0  e  Rn  in  the  linearized 
problem  is  equivalent  to  eliminating  n  of  the  2n+l  rows  of  the  Primal 
Problem  constraints  (22).  All  the  techniques  developed  for  the  Dual  Prob¬ 
lem  (23)  -  (2A)  simplify  slightly  when  applied  to  the  dual  of  this  modified 
problem.  Consequently  the  modified  dual  problem  is  smaller,  easier  to 
solve,  and  requires  less  storage  than  the  unmodified  dual. 
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